function out = simulate_X( pars, Nu, InitX, Nt, localsimM, Phi, Xvec, idx )

[a,thresh]                      = deal(NaN(Nt,Nu.Nsim));
X                               = deal(NaN(Nt+1,Nu.Nsim));
X(1:idx,:)                      = repmat(InitX',[idx 1]);

for nn=idx:Nt
        thresh(nn,:)                      = interp1(Xvec,Phi,X(nn,:));
        a(nn,:)                           = ( localsimM(nn,:) <= thresh(nn,:) );
        X(nn+1,:)                         = (1-pars.k*Nu.Deltat)*X(nn,:) + pars.k*Nu.Deltat*a(nn,:);
end

out.a                           = a;
out.X                           = X;
out.thresh                      = thresh;
    
end